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Abstract 

In this paper, a standard PDE for the pricing of arithmetic average strike Asian call option is pre- 
\ sented. A Crank-Nicolson Implicit Method and a Higher Order Compact finite difference scheme for 

£H | this pricing problem is derived. Both these schemes were implemented for various values of risk free 

rate and volatility. The option prices for the same set of values of risk free rate and volatility was also 
computed using Monte Carlo simulation. The comparative results of the two numerical PDE methods 
shows close match with the Monte Carlo results, with the Higher Order Compact scheme exhibiting a 
better match. To the best of our knowledge, this is the first work to use the numerical PDE approach for 
' pricing Asian call options with average strike. 
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1 Introduction 

Options are one of the most common financial derivatives that are traded both in exchanges as well as over 

>- ■ 

Q\ ' the counter CD 121 [3j HI- The two most common options are the European and the American options. The 
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purchase or sale of the underlying asset for the option takes place at the discretion of the holder of the option 
for a price called the exercise price or the strike price at a fixed time called expiration time (denoted by T) in 
case of European option and any time in [0, T] in case of American option. Options can be further classified 
as call or put depending whether the holder has the right to buy or the right to sell the underlying asset. In 
addition there are numerous other options that can broadly be classified as exotic options flU El El El - The 
specialty of these kinds of options is that the final payoff is more sophisticated and sometimes depends on 
some function of the path of prices of underlying asset. One of the common path dependent exotic options 
is the Asian option H]|2l[5]|6l|71, wnose P av °ff depends on the average historical stock prices. In this paper, 
we will focus on pricing of an Asian call option with arithmetic average strike. 

We begin with the assumption that the stock prices S(t) follow a geometric Brownian motion given by 
the stochastic differential equation HUH [2 EH, 

dS(t) = rS{t)dt + aS{t)dW(t), 

where r is the average rate of growth of asset prices or drift, a is the volatility and W(t)(0 < t < T) is a 
Brownian motion or Wiener process under the risk neutral measure P. The payoff for an Asian call option 
IU|2l|5]|6j|71, with arithmetic average strike is given by, 

V(T) = max ^S{T) - ^ J S(u)du, 



The price of the option at time t G [0, T] (with filtration Ft) is given by the risk-neutral pricing formula 1171 , 

V{t) = E (e- r ^ T -^V{T)\F t 



In order to deal with the challenge of the payoff V(T) being path-dependent we introduce a stochastic 
process I(t) [HEl given by, 

t 

I(t) = [ S(u)du. 





The stochastic differential equation for evolution of I(t) is thus given by, 

dl(t) = S(t)dt. 

Thus the Asian call option price V(S, I, t) for continuous arithmetic average strike satisfies the backward 

pdeehei, 
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-m +rS ds + r S ds^ + s m- rV = - 

Note that the above problem is three dimensional which leads to greater computational costs. This moti- 
vates the reduction of the problem into lower dimension. Q. For this purpose, a new variable R(t) = 
5^tj Jo S( u )du is defined Gl[6]]. This in turn motivates the ansatz V(S, I,t) = S ■ H(R, t) for some func- 
tion H(R, t). It can be shown that the SDE satisfied by R[t) is O, 

dR(t) = (1 + [<J 2 - n)R(t)) dt - aR(t)dW{t). 

Also the function H(R, t) satisfies the PDE 0, 

— + -a 2 R 2 -— + (1 - rR)—- = 0. 

dt 2 dR 2 v ' dR 

The solution of this backward PDE requires a final condition and two boundary conditions which are out- 
lined below HEl. 

i. Final Condition : The final payoff for the option gives the final condition, 

H{R(T),T) = max (1 - ^R(T),0 



ii. Right Hand Boundary Condition : The right hand boundary condition for R — > oo can be obtained by 
observing that since the integral R(t) is bounded, so S — > for R — > oo. For S — > the option is not 
exercised rendering it's value to be 0. Hence, 

H(R, t) = for R^oo. 



iii. 



Left Hand Boundary Condition : The left hand boundary condition for R — > can be obtained from 
the similarity reduction equation. The term RdH/dR — >■ as R —¥ 0. Assuming that H is bounded it 
follows that the term R 2 d 2 H/dR 2 — > as R — > 0. This leads to the boundary condition, 



The problem to be solved now reduces to, 
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H = , for .R oo. 
Once the solution H(R, t) is obtained the price of the Asian option is determined by, 

V(S(0),R(0),0) =S(0)H(R(0), 0), 

where 5(0) is the initial stock price. For pricing of options (especially exotic options) the standard method 
used is Monte Carlo simulation ||8]|9j[T0l. This involves simulating the paths of the underlying asset and 
calculating the option price based on this path. A large number of such simulations are run and the average 
of the option prices from each simulation is taken to be the option price. Several methodologies have been 
adopted for pricing of options using the numerical PDE approach. Some of the most commonly used ones 
are finite differences of lower order |[TTl[T2l[T3l[T4l and higher order compact schemes for American options 
[15] [161 an d option pricing in stochastic volatility model IfTTl . 

2 Crank Nicolson Implicit Method 

The problem of pricing the average strike Asian call option essentially entails solving for equation (fl]). While 
geometric mean Asian option admits closed form solutions (8), the same is not true in case of arithmetic 
average Asian options. As such one has to seek a solution through numerical methods for PDEs |[T2l . There 
are several articles in literature which dwell upon numerical PDE approach to Asian option pricing. One of 
the first papers to deal with numerical PDE pricing of options is by Rogers and Shi 1 1 Q. In this paper, the 
authors first reduce the problem of solving a parabolic PDE in two variables and present a highly accurate 
lower bound. Zvan et al. lfT3l in their technical report, do an extensive analysis of numerical PDE methods 
of Asian options. They discuss the shortcomings of applying the usual numerical PDE techniques used for 
standard options in case of Asian options. In particular they adapt flux limiting techniques from computa- 
tional fluid dynamics (CFD) to tackle the problem of spurious oscillations that arise in Asian options. Vecer 
Ifl2l provided a numerical implementation of the Asian option pricing problem using the 9 method. Dubois 
and Lelievre lfl4l extend the approach by Rogers and Shi IfTTl and propose a scheme which produced fast 
and accurate results. While all these papers IfTTl [T2l [T3l [T4l do examine the pricing problem from the nu- 
merical PDE point of view, the focus is mostly on fixed strike options. Rogers and Shi IfTTl and Zvan et al. 
|[T3l present some results on average strike put options. 

Our main objective in this paper will be to use Higher Order Compact (HOC) scheme for this purpose. 
We will postpone the discussion on this until the next section. In this section we will present the Crank- 



Nicolson Implicit Method (CNIM) for solving equation CD and compare the results with those obtained by 
Monte Carlo simulation lfT"8l . 

CNIM is obtained by taking the average between Forward-Time Centered-Space method (FTCS) and 
Backward-Time Centered-Space method (BTCS). For this puipose, let us define a finite difference dis- 
cretization of the PDE (equation£l])) with the uniform grid t n = t\ + (n — 1) • At, n = 1 : N + 1 and 
Ri = R\ + (i — 1) ■ Ar, i = 1 : M + 1, where At and Ar are the temporal and spatial mesh size respec- 
tively. The values used in this paper are t\ = and tAr + i = T = 1 (i.e 1 year option) with R\ = and 
Rm+i = 5. Let us also define the variables c(R) = \cr 2 R 2 and d(R) = (1 — rR). H(R,t) at the point 
(Ri,t n ) is denoted by Hf. The CNIM discretization of equation CD) is then given by, 
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The above can now be rewritten as, 
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The CNIM can now be written in the matrix form as, 

BH^ = AH {n+1 "> + b {n) 



H M ) T forn = 1 : N + 1. 



(4) 
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From second order finite differences, the one-sided difference is given by, 

—Hi+2 + 4-ffj+i — 3i?i 



dH 
~dR 



2A 



R 



Therefore, applying the same on the left boundary along with backward time approximation we get, 

HI = (1 - 3k)H™ +1 + 4kH% +1 - kHg +1 , 
where k = 2^~- The right boundary H M+l = follows from equation (Q]). The final condition is given by, 



H 



N+l 
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The CNIM formulation above was solved using MatLab ss ™ . The solution was obtained using an iterative 
process which involved a tolerance criterion of |if new (l, 1) — H \d(l, 1)| < e. For the purpose of this 
implementation the number of time and space grid points were taken to be 101 and 501 respectively. The 
tolerance level was taken to be e = 10~ 8 . 



3 Higher Order Compact Scheme 

Finite difference methods have been used for solving ODE and PDE problems for quite a long time. They 
are relatively easier to set up and solve, but require structured mesh |fl9l l20ll . Finite element methods on 
the contrary are more sophisticated, works well with irregular domains and are amenable to unstructured 
meshes. Finite element methods are relatively more challenging to implement. Standard finite difference 
schemes like the CNIM (used in the previous section) are second order accurate. However, in financial 
applications like option pricing a higher level of accuracy is desirable. A direct extension of the central 
difference schemes to achieve higher order accuracy would involve more node points on the stencil. An 
innovative way of achieving higher accuracy with lesser number of nodes in the stencil came by the way of 
higher order compact (HOC) schemes. Spotz and Carey |[T9l and Spotz EDI provide an excellent discourse 
on application of this scheme in case of viscous flow and computational mechanics. 
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0.05 


3.5025 
(0.375339) 


3.1509 
(4.511535) 


5.1148 
(0.365693) 


4.8734 
(4.492267) 


9.3988 
(0.344431) 


9.3486 
(4.556745) 


0.1 


4.1353 
(0.117604) 


4.0124 
(4.505087) 


5.5629 
(0.11831) 


5.4183 
(4.489719) 


9.5333 
(0.117648) 


9.433 
(4.53052) 


0.2 


6.1337 
(0.117679) 


6.1172 
(4.490006) 


7.2951 
(0.117177) 


7.2625 
(4.481982) 


10.547 
(0.117271) 


10.4894 
(4.530701) 


0.3 


8.3256 
(0.11757) 


8.3155 
(4.490814) 


9.3669 
(0.117839) 


9.3484 
(4.518666) 


12.2035 
(0.118941) 


12.163 
(4.500043) 


0.4 


10.5403 
(0.12389) 


10.5358 
(4.483286) 


11.5081 
(0.117663) 


11.4952 
(4.483686) 


14.0885 
(0.116222) 


14.0581 
(4.490074) 



Table 1: Comparison between Monte Carlo Simulation (MC) and Crank Nicolson Implicit Method (CNIM). 
The values in braces represent the CPU time in seconds. The initial stock price was 5(0) = 100. 



Despite it's enormous potential of application to finance problems, HOC schemes have not been used 
much in this area. Zhao et al. fl2~T1 presented a compact scheme for American option pricing with second 
order accuracy in space. Tangman et al. lfT5l [T6l applied a HOC scheme for the pricing of American put 
option. They do a comparative lfl5l analysis with a non-compact fourth order scheme. In their subsequent 
paper they describe an improvement of a method suggested by Han and Wu |[22l . In this chapter we shall 
apply HOC scheme to the setup (Q]) which is written in the following form |[T8l . 

pp. TT f)JT 

where c(R) and d(R) are as defined in the previous section and g = —^M-- We now define notation 5 
IHJHOlUa as follows, 

df d 2 f 

6Rf:= dR> 52Rf:= dm andsoon - 

Thus rewriting Equation (f5]) in terms of finite difference discretization we get, 

Ci5 2 R Hi + di5 R Hi = gi (6) 

In the HOC scheme we derive the leading truncation error terms in terms of finite difference equation making 
use of the original equation. Denoting these terms by r, the HOC scheme is obtained by subtracting r back 
to the original finite difference discretization. Thus we have, 

Ci5 2 R Hi + di5 R Hi -Ti = gi (7) 



where, 
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From the initial PDE £[]) we get, 



dc d 2 H d 3 H dd dH d 2 H 
+ c—-T + —— + d- 



dRdR 2 dR 3 OR OR OR 2 OR 
Therefore this can be rewritten in 5 notation as, 



d 3 H 



OR 3 



1 



5 R Ci5 2 R Hi + 5 R di5 R Hi + di5 2 R Hi 



+ — $R9i 

Ci 



Differentiating the PDE © w.r.t. R once more and simplifying we get, 
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Making use of approximations in equation (fTOl and equation (fT2~l) in the truncation error (n) (equation (HJ) 
and hence substituting t- l in equation (O we obtain, 
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Note that, q = R 2 =>• 5 R Ci = a Ri, 5 R Ci = a 2 and di = 1 — ri?; =>• = -r, <5jj,(fj = 0. 
We define, 
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We now apply the HOC scheme to the above equation (recalling that g = —Qjf) and obtain, 
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The above can now be rewritten as, 



G t H? +1 + KiH? + JiHIL! = A^r+Y + EiH? +L + Fi?f+ 



rn+l 



rn+1 



(15) 



where, 
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The implicit method can be written in the matrix form, 

BH {n) = AH {n+l) + b {n \ 

where H^ n \ A, B and has already been defined in the previous section. As with the case of CNIM the 
number of time and space grid points were taken to be 101 and 501 respectively along with the tolerance 
level of e = 10 -8 . The scheme was implemented using MatLab ss ™ . 
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(0.127529) 


4.0124 
(4.505087) 
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(0.130137) 


5.4183 
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9.4385 
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(4.53052) 


0.2 


5.9919 
(0.121434) 


6.1172 
(4.490006) 


7.1641 
(0.121898) 


7.2625 
(4.481982) 


10.4486 
(0.125619) 


10.4894 
(4.530701) 


0.3 


8.2462 
(0.123275) 


8.3155 
(4.490814) 


9.2902 
(0.125051) 


9.3484 
(4.518666) 


12.1361 
(0.125488) 


12.163 
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0.4 


10.4921 
(0.12519) 


10.5358 
(4.483286) 


11.4607 
(0.120896) 


11.4952 
(4.483686) 


14.0444 
(0.120095) 


14.0581 
(4.490074) 



Table 2: Comparison between Monte Carlo Simulation (MC) and Higher Order Compact Scheme (HOC). 
The values in braces represent the CPU time in seconds. The initial stock price was 5(0) = 100. 



4 Results and Discussion 

In this section we discuss the results obtained by using the CNIM and the HOC schemes as outlined in the 
previous two sections. As already noted we could not find any results for average strike Asian call option 
using numerical PDE methods. For the purpose of comparison we used the Monte Carlo (MC) simulation 
as the benchmark value. We generated the path of a stock prices using the geometric Brownian motion 



process. We generated 50000 such paths and determined the option price from each of the paths generated. 
The average of all these option prices was taken to be the option price, for the purpose of comparison 
with the PDE methods. We generated the option prices using all the three methods for three values of 
r = 0.06, 0.1, 0.2 and five values of a = 0.05, 0.1, 0.2, 0.3, 0.4. 

A comparative study of results from the CNIM and the MC methods showed a close match. The com- 
parative results are presented in Table Q] along with the CPU time in seconds. The option prices for the three 
values of r against the five values of a have been presented in the graphical form in Figures (Q]), (f2]) and ((3]). 
For r = 0.06 (Figure £[])), the match was very close except the case where a = 0.05. This slight difference 
in the option price is reduced when r = 0.1 (Figure ©). The other values for r = 0.1 showed a close match. 
The results were similar for the case r = 0.2 except for a very minimal difference in the case of a = 0.1 
(Figure ([3])). The CPU time in case of CNIM was however considerably lower (< 0.5 seconds) as compared 
with the Monte Carlo simulation (> 4 seconds). 

The results and comparison of the HOC scheme and the MC method indicates excellent consonance. 
A comparison of the results from these two methods in terms of values and CPU time in seconds have 
been presented in a tabular form in Table [2] The option prices from both the methods are very close to 
each other. In fact the results obtained from the HOC scheme show a better match with the MC simulation 
results as compared with the CNIM method. This holds for all the values of r and a and is evident from 
the comparative figures (Figures ©, ©, ©) of HOC and MC. As was the case with CNIM, the CPU time 
taken in case of the HOC scheme is significantly less (< 0.5 seconds) in contrast to Monte Carlo simulation 
which requires at least 4 seconds. 

5 Conclusion 

In this paper we examined several ways of computing the price of an average strike Asian call option, 
namely Monte Carlo simulation and the numerical PDE approach. In case of option pricing, the benchmark 
generally used is Monte Carlo simulation which suffers from some severe drawbacks like computational 
costs and a certain amount of uncertainty of pricing. In contrast, the usage of numerical PDE approaches 
that we have taken results in lesser computational costs and also provides an unique answer. The numerical 
PDE approach in pricing the average strike Asian call option is by and large an unexplored area, since this 
approach applied to Asian option is mostly concentrated on the case of fixed strike. In this paper, we take 
the PDE approach to the pricing problem and present two schemes to accomplish this numerically. Firstly, 
we use the Crank-Nicolson Implicit Method (which is second order) to solve the PDE and hence price the 
option. Then, we present a Higher Order Compact scheme (fourth order) to solve the problem. Finally 
we make a comparison of results obtained from the PDE approach with that of Monte Carlo. The results 
obtained using the two PDE techniques were in excellent agreement with the Monte Carlo results. The 
results obtained using Higher Order Scheme are closer to the Monte Carlo results as opposed to Crank- 
Nicolson Implicit method vis-a-vis Monte Carlo. This is more so in case of lower values of a. To the best of 
our knowledge, this is the first work to use the numerical PDE approach for pricing Asian call options with 
average strike. We believe this work would find more applications in the area of option pricing through the 
PDE approach. 
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Figure 1: Comparison between Monte Carlo Simulation (MC) and Crank Nicolson Implicit Method (CNIM) 
for r=0.06 
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Figure 2: Comparison between Monte Carlo Simulation (MC) and Crank Nicolson Implicit Method (CNIM) 
forr=0.1 



Comparison between Crank Nicolson and Monte Carlo for r=0.2 
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Figure 3: Comparison between Monte Carlo Simulation (MC) and Crank Nicolson Implicit Method (CNIM) 
for r=0.2 



Comparison between Higher Order Compact and Monte Carlo for r=0.06 
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Figure 4: Comparison between Monte Carlo Simulation (MC) and Higher Order Compact Scheme (HOC) 
for r=0.06 



Comparison between Higher Order Compact and Monte Carlo for r=0.1 
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Figure 5: Comparison between Monte Carlo Simulation (MC) and Higher Order Compact Scheme (HOC) 
forr=0.1 



Comparison between Higher Order Compact and Monte Carlo for r=0.2 
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Figure 6: Comparison between Monte Carlo Simulation (MC) and Higher Order Compact Scheme (HOC) 
for r=0.2 



